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It is shown, under weak conditions, that the dynamical evolution of an important class of large 
systems of globally coupled, heterogeneous frequency, phase oscillators is, in an appropriate physical 
sense, time-asymptotically attracted toward a reduced manifold of system states. This manifold, 
which is invariant under the system evolution, was previously known and used to facilitate the 
discovery of attractors and bifurcations of such systems. The result of this paper establishes that 
attractors for the order parameter dynamics obtained by restriction to this reduced manifold are, 
in fact, the only such attractors of the full system. Thus all long time dynamical behavior of the 
order parameters of these systems can be obtained by restriction to the reduced manifold. 
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Systems consisting of many coupled phase os- 
cillators have been used to model diverse situa- 
tions ranging from Josephson junction circuits, to 
circadian rhythms, to synchronization of cardiac 
pacemaker cells. In previous work by us [151 ], it 
was shown that a large class of such models pos- 
sess solutions on an invariant manifold M. It has 
since proved possible to simply obtain various at- 
tractors of the dynamics on M. A remaining open 
question is that of whether such attractors for the 
dynamics on M are also attractors for the dynam- 
ics of the full system, and, if so, whether all of 
the attractors of the full system lie on M. In this 
paper we prove, under very general conditions, 
that, in an appropriate sense, the answer to these 
questions is yes. This result establishes that re- 
striction of consideration to the manifold M can 
be used as an effective computational and analy- 
sis method for obtaining all the typical, long-time 
dynamical behavior of these systems. 



I. INTRODUCTION 

Large systems of coupled phase oscillators with het- 
erogeneous frequency distributions are of general inter- 
est and are the essential modeling tool in past analyses 
of a variety of interesting situations in physics, chem- 
istry, biology, etc. Perhaps the simplest and best known 
such system is the Kuramoto model [l| , which treats the 
synchronization of globally (all-to-all) coupled phase os- 
cillators for which the coupling between pairs of oscilla- 
tors appears as the sine of the phase difference between 
the oscillators. Examples where this basic framework 
has been extended to more complex situations include 
Josephson junction circuits Q, pedestrian induced oscil- 
lation of footbridges 0, E| , systems with time-dependent 
coupling [f| , driven systems describing circadian rhythm 
in mammals Ja, l7|, the effect of time-delay in oscillator 
interactions [g, |9j, the effect of non-unimodal distribu- 
tion of the natural frequencies of the phase oscillators 
[Tol . [Tl| . "communities" of phase oscillators interacting 



with multiple other phase oscillator communities fl2l . , 
the "chimera" model of certain mammals that experience 
sleep with only one of their two brain hemispheres at a 
time [II 03], etc. 

The large number of interesting applications of phase 
oscillator models motivates the attempt to find general 
analysis tools applicable to these models. In this vein, 
it has recently been shown that, in the continuum limit 
(i.e., the number of oscillators approaches infinity), such 
models possess solutions on a reduced manifold of sys- 
tem states [l"> . Furthermore, for the case of a Lorentzian 
distribution of oscillator frequencies, the dynamics on the 
reduced manifold is typically describable by a finite num- 
ber of ordinary differential equations. This finding has 
been utilized to determine attractors and their bifurca- 
tions on the reduced manifold for all the app lications 
previously mentioned (see Refs. @, 1 B 0, H El El El ) • 
Two basic questions remain: (i) are attractors for the 
dynamics restricted to the reduced manifold also attrac- 
tors of the full system; and (ii) are there attractors of 
the full system that do not lie on the reduced manifold? 
Indications of results so far are mixed. On the one hand, 
numerical results from Refs.[y,l9|, and especially [rjj, are 
consistent with the supposition that all attractors of the 
full system lie on the reduced manifold. On the other 
hand, Refs. [HI, find long-time asymptotic behavior 
that is not on the reduced manifold. The result of our 
paper is that, in an appropriate sense (that we specify 
later in this paper), the reduced manifold is globally at- 
tracting provided that the spread A in the distribution of 
oscillator frequencies is nonzero. In particular, for A > 
all attractors of the full system lie on the reduced man- 
ifold, and all attractors of the dynamics on the reduced 
manifold are attractors of the full system. This greatly 
facilitates the task of finding the attractors of the full 
system, since they now can be sought using the reduced 
system. The result also resolves the puzzle posed by the 
previous results, since the finding by Refs. [£3, 03] 

of long 

time motion not on the reduced manifold was for the case 
of A = 0, while the opposite indication from the numer- 
ical results of Refs. [J ||J HH treated situations in which 
A > 0. 
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II. FORMULATION 

We begin by noting that the models in the class of 
problems in which we are interested all involve the deter- 
mination of a distribution function F '(0,uj,t), where 9 is 
the phase of an oscillator, and ui is the natural frequency 
an oscillator would have in isolation from the outside 
world (e.g., from other oscillators); Fdddui is the fraction 
of oscillators at time t whose phases and natural frequen- 
cies lie in the range [9, 9 + d9] and [u>, lu + du>]. Since the 
natural frequency uj of an oscillator is assumed not to 
change with time, the marginal frequency distribution, 



F(6,w,t)d0, 



(1) 



is time independent. The key quantity characterizing 
the macroscopic behavior of the distribution function F 
is the "order parameter" r(t) originally introduced by 
Kuramoto [1| and defined by 



r(t) 
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F(9,uj,t)e- ie d6duj. 



(2) 



Since the number of oscillators is conserved, F obeys an 
oscillator continuity equation, 



■at + W {veF) 



o. 



(3) 



and, for all of the problems previously mentioned 
(Refs.[|[ to [3]), vg(9,t) is expressible in the form 0, 



several distribution functions, i.e., F and H in ([3]) and ((4]) 
are replaced by F G and H a (a — 1, 2, • • • , s, where s is the 
number of distributions) and each H a is a function of all 
the order parameters r-y , r% , • • • , r s ; in the case of pedes- 
trian induced oscillation of footbridges 0,13] > H = ky(t), 
where y(t) is the side-to-side acceleration of the bridge, 
which obeys a damped oscillator equation driven by r(t), 
where r(t) represents the effects of the pedestrians. 

For a general problem of the type described above, 
as the system evolves, H(t) will change self consistently 
with t. We will show in what follows that, whatever is 
the evolution of H(t), in the long time limit, solutions 
for the order parameter evolution r(t) (or r a (t)) obey 
the differential equation that applies for evolution on the 
reduced manifold of Ref.[l5|. Because the precise time 
dependence of H (t) will not matter in the derivation of 
our result, it suffices to consider H(t) as some general 
function of time without regard to how this time depen- 
dence is determined. 

Expanding the distribution F as a Fourier series in 9, 
we write F in the form, 



F(9,u,t) 



2tt 



[l+F + (9,u,t) + F_(6,uj,t)], (5) 



F+(9,u,t) = J2F n (uj,t)e 

n=l 

oo 



(6) 



v e (0,u,t) =^ + - [H{t)e- M -H*(t)e ie ] . (4) 

Equations ([3]) and ((H) constitute an w-dependent par- 
tial differential equation in the two real variables (9, t) 
to be solved subject to an initial condition F(6,u>,0). 
The problem is ostensibly complicated by the fact that 
the time dependence of the quantity H may, in general 
(Refs.[l| to (HI), depend on F through the complex order 
parameter r(t) defined by @, as well as through other, 
non-phase-oscillator variables obeying auxiliary dynam- 
ical equations, which themselves may @, Q depend on 
r(t), or (asinRef. @) through explicit external time de- 
pendence of system parameters. Here are some examples: 
for the classical Kuramoto [l| problem, H = kr(t), where 
k is the strength of the coupling between oscillators; for 
the circadian rhythm problem H = kr(t)+A, where 
A is a constant reflecting the strength of the diurnal drive 
of the day-night sunlight cycle, and A might be given an 
explicit time dependence, A = A(t), to model variation 
between sunny and cloudy days; for the case of time de- 
lay in the response of oscillators to other oscillators in the 
system H = k p(r)r(t — r)dr, where p(r) is the 

distribution function [9( of delays along the links between 
oscillators; in the cases treated in Refs.[l2l. [l3j (commu- 
nities of oscillator s'). Illl (nonunimodal frequency distri- 
bution g(u))), and [lj;[lJ| (the chimera model), there are 



We note that J Q 27r F±d9 = and that, assuming absolute 
convergence of the Fourier series, the analytic continua- 
tion of F + (F_) into Im{9) > (Im(9) < 0) has no singu- 
larities and decays exponentially to zero as Im{9) — > +oo 
(Im(9) — > — oo). As will soon become evident, the de- 
composition of F given by ([5]) is a key step. We note 
that since F_ = F+ for real (9,ui), it suffices to consider 
only F + . Substituting ([5]) into Eqs. (|3l4[ ) and projecting 
the result onto the function space spanned by the basis 
functions, e %e , e 2l ° , e ZlB , • • • , we obtain 



dF+ 



d_ 

89 
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H" 
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F, 
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As previously noted, our result will not depend on the 
precise time dependence of H(t). Thus, whatever is the 
dependence of H(t), we can formally regard it as given. 
Adopting this viewpoint, Eq.((7]) is linear in F+ with 
an inhomogeneous driving term on the right hand side 
(namely, \H* exp(i6*)). As such, we can write F + as 



F, 



F,+F, 



(8) 



where F+ is a homogeneous solution to {7J and F + is an 
inhomogeneous solution. An inhomogeneous solution is 
given by taking the Fourier coefficients of F+ to be given 



3 



by F n (u>,t) = [a(uj, t)] n , as proposed in Ref . [3~5I] . When 
this ansatz is used in ([7]), it is found that ([7]) is indeed 
satisfied if a (a;, i) satisfies 



da 
~dt 



(Ha 2 -H*) = 0, 



(9) 



which, for each value of u>, is an ordinary differential 
equation in time t. We note further that, as shown in 

[15| . |a(w,i)| < 1 so that the summation of the Fourier 

- i 

series for F + converges and yields 



what we require is that 



where 



lim f+(9,t)=0, 

t— ► OC 



/+(*,*) 



F + {d,w,t)g{u})du. 



(15) 



(16) 



In what follows we will demonstrate that Eq. (fl5jl indeed 
holds under very general conditions. 



F + = 



1 



(10) 



This form of F+ specifies the reduced manifold found in 
Ref. [l5| . Thus, at long time, F + would tend to F + (i.e., 
F would tend to the reduced manifold) if lim^oo F + — 
0. However, a simple counter examples shows that this 
cannot always be true. In particular, if H = 0, Eq.([7]) 
has homogeneous solutions, 



in{9—ujt) 



(ii) 



for any set {A„} for which this series converges. Since 
the magnitude of each term of the summation in (fTTjl is 
time-independent (u> is real), F + does not go to zero at 
t —> oo. On the other hand, we note that, as t increases, 
the individual terms, e m ( e -" t ) j oscillate more and more 
rapidly in u>. Thus for any such term 



g(cj)A n e in ^-^duj 



decays exponentially in time for sufficiently smooth g(u>). 
For example, our subsequent considerations will be for 
the case of a Lorentzian frequency distribution, 



i 



i 



ttlo 2 + A 2 2i:i\uj-iA ui + iA 



1 



1 



for which 



T — A p 1 



iB—nAt 



, (12) 



(13) 



which decays exponentially to zero as t — > oo provided 
that A > (the case A = corresponds to g(u>) being a 
delta function in ui) . [We remark that the mean value of 
u> has been taken to be zero in (TH?)) , but that no gener- 
ality is lost by this, as a mean value can be restored by 
the change of variables, 9' = 9 — fit, ui' = cu + fL] Thus, 
while we cannot expect to show that F + — > as t — » oo, 
Eq. (fT"3")) suggests that this may not be necessary to obtain 
order parameter dynamics that tend to the order param- 
eter dynamics that applies on the reduced manifold. In 
particular, noting from (O and © that 



r(t) 



2tt 



F + ge- t9 d9duj, 



(14) 



III. DEMONSTRATION OF THE MAIN 
RESULT 

To show that (fT5|) applies, we now assume that the 
analytic continuation of F + (9,uj,t) into Im(u>) < has 
no singularities in Im(w) < and approaches zero as 
Im(u>) — * — oo. To show that this last assumption is a 
consistent one, let \oj\ be very large, \oj\ 3> H . Then the 
homogeneous version of Eq. ([7]) for F+ is approximately 



dF+ 
dt 



dF+ 
' 89 



0. 



which has solutions for its Fourier ^-components F n ~ 
exp[in(9 — cot)] which go to zero as Im(u>) — > — oo. (It 
was to achieve this that we have introduced the decom- 
position of F given by Eqs.([5]) and (J6j) . ) We will further 
discuss this analyticity assumption at the end of this pa- 
per. 

We now specialize to the case of Lorentzian g(u>), 
Eq. p2]) . We multiply the homogeneous version of Eq. 
by g(u>)clio, integrate the result from uj = —R tow = +R, 
analytically continue into the complex w-plane, close the 
integration path with a semicircle of radius R in the lower 
half w-plane, and let R — > oo. Using our assumption that 
F+(9,uj,t) is analytic in the lower half w-plane and de- 
cays to zero as Im(u>) — * — oo, the integral along the large 
semicircle approaches zero as R — ► oo, and the integrals 
from uj = — oo to ui = +oo along the real w-axis may thus 
be evaluated as the residue of the enclosed pole of g(ui) 
at oj = -iA (see Eq.([I2])). This yields 



df+(6,t) , d 



dt 



09 



,t) = -I 



[v(6,t)U(e,t)]=0, (17) 



- ie H(t)-e l9 H*{t)) , (18) 



where f+(9,t) = F+(9,-iA,t). 

We now introduce a conformal transformation of the 
upper half complex 0-planc into the unit disc, z = e lS . 
Equations (fT7|) and (fT5|) then become 



df+(z,t) 
dt 



° [ v{z,t)f+{z,t) 



dz 



0. 



(19) 
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where 

v(z,t) = Az + ^(H(t)- z 2 H*(t)) 
f+(z,t) = f+(9,t)/e ie . 
Noting that (fT9|) can be written as 



0. 



df+(z,t) ~ ,,dv{z,t) 



dt 



dz 



(20) 
(21) 

(22) 



where d/dt = d/dt + vd/dz, we can integrate (|2"2l along 
the characteristics of this equation to obtain 



/+(*, t) = f+(Z(z, 0), 0) exp i-r,(z, t)] , (23) 



where 



r)(z,t) = 



dv(z',t') 
dz' 



dt', 



z'=Z(z,t') 



and the characteristics are given by the orbit equation 
dZ{z, t') 



dt' 



= v(Z(z,t'),t'), 



(24) 



with the final condition Z(z,t) = z. Thus Z(z,t') for 
t' < t represents the location Z of the orbit v(z', t') that 
winds up at point z at time t. It is useful to rewrite (124|) 
by introducing 

Z = pe i < t >, H = he ip 

with h, 0, p and <f> real. The real and imaginary parts of 
then give 



P 



dp 

dt 7 

d0 



v p = pA+^(l- i o 2 )cos( 



-(l + p 2 )sin(0-/3). 



P), (25) 
(26) 



We note from (|23|) that when p = 1, we have dp/dt' = 
A > 0. Thus for final conditions on p — 1, the orbits 
backward in time move into p < 1. Thus |Z(z,£')| < 1 
for \z\ < 1 and t' < t (i.e., Z(z,t') is in the unit disc). 
We wish to show that f+(z,t) — > as t — » +oo. From 
(f2"3")l we see that this will be the case if 

lim Re[r)(z,t)] — +oo. 

t — >oo 

In order to show this, we first note that the real part of 
dv(z' , t')/dz' is simply one half the divergence of the two 
dimensional flow v = v p (p,(f)po + v$(p, 4>)4>o (where v p 
and v<j> are given by (|25|) and (|26|) . and po and (fio are 
unit vectors in the p and <f> directions), i.e., 



Re 



f dv(z',f) 
\ dz~> 



z>=Z(z,t>) 



(27) 



Equation ([2T|) is most easily demonstrated in rectangu- 
lar co-ordinates: z' = x + iy, v(z',t') — v x (x,y,t') + 



iv y (x,y,t'), where v x ,v y ,x and y are real. Then (f27|) 
immediately follows by setting v = v x x + v y yo, and us- 
ing the Cauchy-Riemann condition, dv x /dx — dv y /dy, 
in the expression for the divergence in rectangular coor- 
dinates, V • v = dv x /dx + dvy/dy. Now evaluating V • v 
in polar coordinates (p,4>)i we have 



Vv 



ld_ 

pdp 



1 dv,p 



2[A-hpcos((f> -13)}. (28) 



Solving (|28p for hcos(4>—P) in terms of V-v and inserting 
the result in Eq. (P5|) for dp/dt', we obtain after some 
rearrangement 



Re 



( dv{z',t') 
\ dz' 



J + P 2 (z,t') 

z'=Z(z,V) l — p 2 {z,t') 



=A 



+ — Hl-p\z,t')}. (29) 

Inserting ([2^)1 into the integral for r/(z,t) and choosing a 
fixed reference time T satisfying < T < t, we have 



Re[r)(z,t)] 



i - i V dz 

2 



dv(z',t') 



In 

■A 



Re 

l-p 2 (z,t-T ) 
l-p 2 (z,0) 

„ l- P 2 (z,t') 



dt' 



z'=Z(z,f) 



dt'. 



(30) 



We are interested in final (t' = t) conditions on the unit 
circle, Z(z, t) = z = e l8 for 8 real, and their continuation 
into the unit disc \z\ < 1, corresponding to p < 1 at the 
final time t' = t. For p sufficiently near one, Eq. (|2"5l) 
shows that dp/dt' = A. Thus by the continuity of the 
right hand side of (|25|) . there is an annulus in the Z-plane, 
1 > P ^ Po> m which dp/dt' > 0, implying that as t' is 
reduced from t (i.e., t — t' is increased), p moves uniformly 
from p = 1 at time t to smaller values. Thus any final 
point in the annulus eventually enters the disc p < po < 1 
and never leaves it. We can therefore choose the time T 
such that for all final conditions \Z(z,t)\ = \z\ < 1, we 
have 



p{z,t')<p(z,t-T)<p <l 



(31) 



where the first inequality applies for < t' < t — T. We 
consider T to be held fixed, and we ask how r](z, t) be- 
haves as t — > +oo. By (f2"8")) the integrand in the first 
of the three terms of (|30|) is bounded, and, since the in- 
tegration range in t' (namely, T) for this term is fixed, 
we conclude that the first term is bounded. By (|3"l"T) the 
second term in ([3H)l is also bounded. Again by (|3"Tj) the 
integrand of the third term of (f3T)|) is positive and greater 
than one. Thus the third term exceeds (t — T)A. Hence 
Re[rj(z,t)] — > +oo for t — > +oo, if A > 0, thus demon- 
strating that the exponential factor in (|23p goes to zero 
for large time [l7| . We therefore conclude from (f3T)]) and 



5 



(123)) that Eq.lJT&J) is satisfied if A > 0, which is the desired 
result. 

In particular, at large t the order parameter will ap- 
proach the quantity f_ a(u),t)g(u))du, where a(u>,t) 
evolves by Eq.Q. This implies that r(t) will satisfy the 
differential equation, 

^ + Ar(t) + \ [H(t)r\t) H*(t)} = 0, (32) 

which follows from multiplying ^ by g(u))du), integrating 
from u> = — oo to uj = +oo and, as done previously, us- 
ing the residue method to evaluate the integrals. Hence, 
for A > 0, the long time dynamics of the order parame- 
ter r(t) is governed by the ordinary differential equation 
fEq. ([3"2")) and Ref. [HI) that describes its dynamics for dis- 
tribution functions F on the reduced manifold. This is 
our main result. 



IV. DISCUSSION 

The principal conditions for the applicability of our 
result are that the initial condition is such that, when 
F + {9, u>, t) is continued into the complex w-plane, it is an- 
alytic in Im(uj) < and decays to zero as Jm(w) — ► — oo. 
As discussed in Ref. 15], if these conditions are satisfied 
initially, then they are also satisfied for all t > 0. What 
happens if the condition at Im{uj) — > — oo is not satisfied 
initially? Here, a simple example [l8| may be instructive. 
We again consider the case H = 0. Say the initial condi- 
tion on F + has a component exp(m# + i^/Lu) with 7 real 
and positive. This initial condition violates our assump- 
tion of decay to zero as Im(u>) — > —00. However, use of 
this initial condition in Eqs. (|3l4p with H = yields the 
solution exp(in6> + 2(7 — t)u>), and, for large enough time, 
t > 7, the result satisfies the required condition that it 
approaches zero as Im(u>) — > —00 [19j. Thus, even if 



our desired condition at Im(uj) — > — 00 is not satisfied 
initially, in many cases, the result that the long time dy- 
namics of r(t) is described by Eq. (f3"2"|) may still apply. 

We now connect our result with the concept of an in- 
ertial manifold. An inertial manifold M with respect to a 
distance metric /i satisfies the condition that, for any ini- 
tial condition in the state space, the subsequent system 
evolution is such that the distance between the evolved 
orbit and the manifold M as measured by the metric jj, 
approaches zero as t — > +00. What we have shown in 
this paper is that, in the space of distribution functions 
F{9, u>, t), our reduced manifold (Eq. pH)) ) is inertial with 
respect to the proper distance metric fi. In particular, 
this is so if we take the distance between two distribu- 
tion functions Fi(9,ui,t) and F2(0,u>,t) to be defined by 



fJ,(F 1 ,F 2 




- F 2 )g(uj)duj 



1/2 



dO 



(33) 

(For Fi not on the reduced manifold M, the distance 
from Fi to M is n{F\,F2) minimized over all F 2 on M .) 
Note that, by this choice of distance metric, the problem 
associated with the example of Eq. fTTj) is avoided. 

Finally, we remark that, while our result is for the 
special case of a Lorentzian distribution of oscillator fre- 
quencies fEq. (fT2"]0 . we believe that this restriction does 
not greatly limit the usefulness of the resulting formu- 
lations for discovering typical svstem behavior. Indeed, 
past numerical experiments 0, fill ] comparing results ob- 
tained using Lorentzian g(u>) and using Gaussian g(uj) 
were found to yield qualitatively identical bifurcation 
structures. 
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